knitr::opts_chunk$set(warning=FALSE, message=FALSE)
options(scipen=0, digits=4)Below are code and outputs for the ROI-based analyses in the submitted manuscript: “Neural substrates for moral judgments of psychological versus physical harm”.
To see larger versions of any figure, right-click, copy image location, and paste the address to a new tab on your browser.
If you have any questions and/or comments, please email Lily Tsoi: lily [dot] tsoi [at] bc [dot] edu.
Install packages and load libraries
packages <- c("rmarkdown", "knitr", "tidyverse", "broom", "lme4", "ordinal", "lsmeans")
packages_new <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(packages_new)) install.packages(packages_new)
lapply(packages,library,character.only=T)Data files can be found on GitHub: https://github.com/tsoices/psych-phys-harm
Analyses require the following files:
Make sure these files are in the same directory.
files <- c("ROI_PSCs.csv", "ROI_MVPA.csv")
dat_names <- c("dat_psc_orig", "dat_mvpa_orig")
for(i in 1:length(files)) {
assign(dat_names[i], read.csv(paste(params$directory, files[i], sep='/')))
}
# change the order of levels
dat_psc_orig$Group <- factor(dat_psc_orig$Group, levels=c("NT", "ASD"))
dat_psc_orig$Violation <- factor(dat_psc_orig$Violation, levels=c("PH", "PS", "N"))Analyses are based on the following:
Examining behavioral responses in the scanner
Data organization
Organize behavioral data
# calculate mean rating as variable on y-axis
dat_behav <- dat_psc_orig%>%
filter(Violation == 'PH' | Violation == 'PS' | Violation == 'N') %>%
group_by(Subject, Violation, Group, Item, Key) %>%
summarise(mean=mean(Key)) %>%
droplevels.data.frame(.)
dat_behav$Item <- match(dat_behav$Item, unique(sort(dat_behav$Item))) # ordering items such that it doesn't care about purity itemsggplot(dat_behav, aes(x=Violation, y=mean, color=Violation)) +
stat_summary(fun.data="mean_cl_boot", position=position_dodge(0.2), size=1) +
ylim(1,4) +
facet_wrap(~Group, ncol=2, labeller=labeller(Group=c(NT="Neurotypical", ASD="ASD")), scales="free_y") +
scale_x_discrete(labels=c('Physical','Psychological', 'Neutral')) +
scale_colour_manual(name="Condition", labels=c("Physical", "Psychological", "Neutral"), values=c("red", "darkorchid4", "slategray")) +
ylab("Rating\n(1=not at all, 4=very)") +
xlab("Condition") +
theme_bw() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=16,face="bold"),
axis.title.y=element_text(margin=margin(r=18)),
axis.title.x=element_text(margin=margin(t=18)),
plot.title=element_text(size=18,face="bold", margin=margin(b=20), hjust=0.5),
legend.text=element_text(size=14),
legend.title=element_text(size=14,face="bold"),
strip.text=element_text(size=14),
panel.grid.major.x=element_blank(),
panel.grid.major.y=element_blank())ggplot(dat_behav, aes(y=mean, x=Item, color=Violation)) +
stat_summary(fun.data="mean_cl_boot", na.rm=TRUE) +
ylim(1,4) +
facet_wrap(~Group, ncol=2, labeller=labeller(Group=c(NT="Neurotypical", ASD="ASD")), scales="free_y") +
scale_colour_manual(name="Condition", labels=c("Physical", "Psychological", "Neutral"), values=c("red", "darkorchid4", "slategray")) +
ylab("Rating\n(1=not at all, 4=very)") +
xlab("Item") +
theme_bw() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=16,face="bold"),
axis.title.y=element_text(margin=margin(r=18)),
axis.title.x=element_text(margin=margin(t=18)),
plot.title=element_text(size=18,face="bold", margin=margin(b=20), hjust=0.5),
legend.text=element_text(size=14),
legend.title=element_text(size=14,face="bold"),
strip.text=element_text(size=14),
panel.grid.major.x=element_blank(),
panel.grid.major.y=element_blank())Define the model
dat_behav$Key <- factor(dat_behav$Key)
model_behav <- clmm(Key ~ Violation*Group + (1|Subject) + (1|Item),
data=dat_behav,
link="probit",
na.action=na.omit)Test interaction between Group and Condition
anova(model_behav, update(model_behav, . ~ . -Violation:Group))Likelihood ratio tests of cumulative link models:
formula: link: threshold:
update(model_behav, . ~ . - Violation:Group) Key ~ Violation + Group + (1 | Subject) + (1 | Item) probit flexible
model_behav Key ~ Violation * Group + (1 | Subject) + (1 | Item) probit flexible
no.par AIC logLik LR.stat df Pr(>Chisq)
update(model_behav, . ~ . - Violation:Group) 9 2443 -1213
model_behav 11 2445 -1212 2.07 2 0.36
Test main effect of condition
anova(update(model_behav, . ~ . - Violation:Group), update(model_behav, . ~ . - Violation:Group - Violation))Likelihood ratio tests of cumulative link models:
formula: link: threshold:
update(model_behav, . ~ . - Violation:Group - Violation) Key ~ Group + (1 | Subject) + (1 | Item) probit flexible
update(model_behav, . ~ . - Violation:Group) Key ~ Violation + Group + (1 | Subject) + (1 | Item) probit flexible
no.par AIC logLik LR.stat df Pr(>Chisq)
update(model_behav, . ~ . - Violation:Group - Violation) 7 2515 -1251
update(model_behav, . ~ . - Violation:Group) 9 2443 -1213 76 2 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
lsmeans(model_behav, pairwise ~ Violation)NOTE: Results may be misleading due to involvement in interactions
$lsmeans
Violation lsmean SE df asymp.LCL asymp.UCL
PH 1.0778 0.1362 NA 0.8108 1.3448
PS 0.7947 0.1350 NA 0.5300 1.0594
N -0.9417 0.1421 NA -1.2203 -0.6631
Results are averaged over the levels of: Group
Confidence level used: 0.95
$contrasts
contrast estimate SE df z.ratio p.value
PH - PS 0.2831 0.1347 NA 2.102 0.0895
PH - N 2.0195 0.1534 NA 13.163 <.0001
PS - N 1.7364 0.1519 NA 11.432 <.0001
Results are averaged over the levels of: Group
P value adjustment: tukey method for comparing a family of 3 estimates
Data organization
Analyses described here are over the entire time course. Future versions of this document will let you select a time window and automatically refresh the outputs related to that time window.
Time points of interest (seconds):
time_entire <- dat_psc_orig %>% filter(Timepoint >= 6 & Timepoint <= 26)
time_background <- dat_psc_orig %>% filter(Timepoint >= 6 & Timepoint <= 10)
time_action <- dat_psc_orig %>% filter(Timepoint >= 12 & Timepoint <= 14)
time_outcome <- dat_psc_orig %>% filter(Timepoint >= 16 & Timepoint <= 18)
time_intent <- dat_psc_orig %>% filter(Timepoint >= 20 & Timepoint <= 22)
time_judgment <- dat_psc_orig %>% filter(Timepoint >= 24 & Timepoint <= 26)
dat_psc_long <-
time_entire %>%
filter(Violation == 'PH' | Violation == 'PS' | Violation == 'N') %>%
group_by(Subject, Violation, ROI, Group, Item, Key) %>%
summarise(PSC=mean(PSC)) %>%
droplevels.data.frame(.)
dat_psc_long$Item <- as.factor(dat_psc_long$Item)Subset data to NT group only and define model
data_nt <- subset(dat_psc_long, Group == "NT")
model_nt <- lmer(PSC ~
Violation*ROI +
(1|Subject) + (1|Item), data=data_nt, REML=FALSE)
model_nt_1 <- lmer(PSC ~
Violation*ROI +
(Violation|Subject) + (1|Item), data=data_nt, REML=FALSE)Test the condition x ROI interaction
anova(model_nt, update(model_nt, . ~ . - Violation:ROI))Data: data_nt
Models:
update(model_nt, . ~ . - Violation:ROI): PSC ~ Violation + ROI + (1 | Subject) + (1 | Item)
model_nt: PSC ~ Violation * ROI + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_nt, . ~ . - Violation:ROI) 9 5907 5961 -2945 5889
model_nt 15 5918 6009 -2944 5888 0.96 6 0.99
Test the main effect of condition
anova(update(model_nt, . ~ . - Violation:ROI), update(model_nt, . ~ . - Violation:ROI - Violation))Data: data_nt
Models:
update(model_nt, . ~ . - Violation:ROI - Violation): PSC ~ ROI + (1 | Subject) + (1 | Item)
update(model_nt, . ~ . - Violation:ROI): PSC ~ Violation + ROI + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_nt, . ~ . - Violation:ROI - Violation) 7 5911 5954 -2949 5897
update(model_nt, . ~ . - Violation:ROI) 9 5907 5961 -2945 5889 8.42 2 0.015 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# test pairwise contrasts
lsmeans(model_nt, pairwise ~ Violation)Loading required namespace: lmerTest
NOTE: Results may be misleading due to involvement in interactions
$lsmeans
Violation lsmean SE df lower.CL upper.CL
PH 0.11599 0.03906 52.59 0.03764 0.1943
PS 0.15406 0.03906 52.59 0.07571 0.2324
N 0.02426 0.03906 52.59 -0.05409 0.1026
Results are averaged over the levels of: ROI
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
PH - PS -0.03807 0.04356 36.89 -0.874 0.6600
PH - N 0.09173 0.04356 36.89 2.106 0.1025
PS - N 0.12980 0.04356 36.89 2.979 0.0138
Results are averaged over the levels of: ROI
P value adjustment: tukey method for comparing a family of 3 estimates
Test whether including random slope of condition improves model
anova(model_nt, model_nt_1)Data: data_nt
Models:
model_nt: PSC ~ Violation * ROI + (1 | Subject) + (1 | Item)
model_nt_1: PSC ~ Violation * ROI + (Violation | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
model_nt 15 5918 6009 -2944 5888
model_nt_1 20 5928 6048 -2944 5888 0.55 5 0.99
Figure: time course for each condition by ROI
dat_psc_orig$ROI <- factor(dat_psc_orig$ROI, levels=c("RTPJ", "LTPJ", "PC", "DMPFC"))
dat_psc_orig$Violation <- factor(dat_psc_orig$Violation, levels=c("PS", "PH", "N"))
dat_psc_nt <- dat_psc_orig%>%
filter(Group == 'NT' & (Violation == 'PH' | Violation == 'PS' | Violation == 'N')) %>%
group_by(Subject, Violation, ROI, Timepoint) %>%
summarise(PSC=mean(PSC))
cols <- c("PS"="darkorchid4", "PH"="red", "N"="slategray")
rois <- c(RTPJ="rTPJ", LTPJ="lTPJ", PC="precuneus", DMPFC="dmPFC")
ggplot(dat_psc_nt, aes(y=PSC, x=Timepoint, color=Violation, fill=Violation)) +
geom_smooth(na.rm=TRUE) +
facet_wrap(~ROI, ncol=4, labeller=labeller(ROI=rois)) +
annotate("rect", xmin=5, xmax=27, ymin=-Inf, ymax=Inf, alpha=.1) +
scale_x_continuous(limits=c(0,28), breaks=seq(0,28,2)) +
ylab("Percent signal change (PSC)") +
xlab("Timepoint (s)") +
scale_fill_manual(name="Condition\n", labels=c("Psychological harm", "Physical harm", "Neutral act"), values=cols) +
scale_colour_manual(name="Condition\n", labels=c("Psychological harm", "Physical harm", "Neutral act"), values=cols) +
theme_bw() +
theme(axis.text=element_text(size=16),
axis.title=element_text(size=24,face="bold"),
axis.title.y=element_text(margin=margin(r=20)),
axis.title.x=element_text(margin=margin(t=20)),
legend.text=element_text(size=20),
legend.title=element_text(size=24,face="bold"),
legend.key.size=unit(3, "lines"),
strip.text=element_text(size=28, face="bold"),
panel.grid.major.x=element_blank(),
panel.grid.major.y=element_blank())Subset data by ROI (NT only) and define models
dat_rtpj_nt <- subset(dat_psc_long, ROI == "RTPJ" & Group == "NT")
dat_ltpj_nt <- subset(dat_psc_long, ROI == "LTPJ" & Group == "NT")
dat_pc_nt <- subset(dat_psc_long, ROI == "PC" & Group == "NT")
dat_dmpfc_nt <- subset(dat_psc_long, ROI == "DMPFC" & Group == "NT")
model_rtpj_nt <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
dat=dat_rtpj_nt, REML=FALSE)
model_ltpj_nt <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
dat=dat_ltpj_nt, REML=FALSE)
model_pc_nt <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
dat=dat_pc_nt, REML=FALSE)
model_dmpfc_nt <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
dat=dat_dmpfc_nt, REML=FALSE)rTPJ
# test pairwise contrasts
lsmeans(model_rtpj_nt, pairwise ~ Violation)$lsmeans
Violation lsmean SE df lower.CL upper.CL
PH 0.2007 0.05077 42.85 0.098307 0.3031
PS 0.2216 0.05077 42.85 0.119200 0.3240
N 0.1084 0.05077 42.85 0.006032 0.2108
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
PH - PS -0.02089 0.04785 35.72 -0.437 0.9005
PH - N 0.09228 0.04785 35.72 1.929 0.1454
PS - N 0.11317 0.04785 35.72 2.365 0.0597
P value adjustment: tukey method for comparing a family of 3 estimates
lTPJ
# test pairwise contrasts
lsmeans(model_ltpj_nt, pairwise ~ Violation)$lsmeans
Violation lsmean SE df lower.CL upper.CL
PH 0.2595 0.06081 39.21 0.13651 0.3824
PS 0.2707 0.06081 39.20 0.14768 0.3936
N 0.1401 0.06081 39.21 0.01711 0.2630
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
PH - PS -0.01118 0.07468 35.54 -0.150 0.9877
PH - N 0.11940 0.07468 35.55 1.599 0.2594
PS - N 0.13058 0.07468 35.54 1.748 0.2017
P value adjustment: tukey method for comparing a family of 3 estimates
precuneus
# test pairwise contrasts
lsmeans(model_pc_nt, pairwise ~ Violation)$lsmeans
Violation lsmean SE df lower.CL upper.CL
PH -0.03654 0.04516 49.59 -0.12727 0.05419
PS 0.03257 0.04516 49.59 -0.05817 0.12330
N -0.10693 0.04516 49.59 -0.19767 -0.01620
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
PH - PS -0.0691 0.04934 35.44 -1.401 0.3514
PH - N 0.0704 0.04934 35.44 1.427 0.3382
PS - N 0.1395 0.04934 35.44 2.827 0.0204
P value adjustment: tukey method for comparing a family of 3 estimates
dmPFC
# test pairwise contrasts
lsmeans(model_dmpfc_nt, pairwise ~ Violation)$lsmeans
Violation lsmean SE df lower.CL upper.CL
PH 0.04848 0.07077 33.34 -0.09544 0.1924
PS 0.10088 0.07077 33.34 -0.04304 0.2448
N -0.03575 0.07077 33.34 -0.17967 0.1082
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
PH - PS -0.05240 0.07065 513 -0.742 0.7387
PH - N 0.08423 0.07065 513 1.192 0.4584
PS - N 0.13663 0.07065 513 1.934 0.1302
P value adjustment: tukey method for comparing a family of 3 estimates
Subset data to ASD group only and define model
data_asd <- subset(dat_psc_long, Group == "ASD")
model_asd <- lmer(PSC ~
Violation*ROI +
(1|Subject) + (1|Item), data=data_asd, REML=FALSE)
# model with slope of condition does not converge, so no comparison will be made between model w/ slope and model w/o slope.
# model_asd_1 <- lmer(PSC ~
# Violation*ROI +
# (Violation|Subject) + (1|Item), data=data_asd, REML=FALSE)Test the condition x ROI interaction
anova(model_asd, update(model_asd, . ~ . - Violation:ROI))Data: data_asd
Models:
update(model_asd, . ~ . - Violation:ROI): PSC ~ Violation + ROI + (1 | Subject) + (1 | Item)
model_asd: PSC ~ Violation * ROI + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_asd, . ~ . - Violation:ROI) 9 3506 3556 -1744 3488
model_asd 15 3517 3600 -1744 3487 0.85 6 0.99
Test the main effect of condition
anova(update(model_asd, . ~ . - Violation:ROI), update(model_asd, . ~ . - Violation:ROI - Violation))Data: data_asd
Models:
update(model_asd, . ~ . - Violation:ROI - Violation): PSC ~ ROI + (1 | Subject) + (1 | Item)
update(model_asd, . ~ . - Violation:ROI): PSC ~ Violation + ROI + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_asd, . ~ . - Violation:ROI - Violation) 7 3513 3551 -1749 3499
update(model_asd, . ~ . - Violation:ROI) 9 3506 3556 -1744 3488 10.9 2 0.0042 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# test pairwise contrasts
lsmeans(model_asd, pairwise ~ Violation)NOTE: Results may be misleading due to involvement in interactions
$lsmeans
Violation lsmean SE df lower.CL upper.CL
PH 0.11668 0.03348 37.28 0.04886 0.18451
PS 0.17087 0.03348 37.28 0.10304 0.23870
N 0.02851 0.03348 37.28 -0.03932 0.09633
Results are averaged over the levels of: ROI
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
PH - PS -0.05419 0.04062 39.6 -1.334 0.3852
PH - N 0.08818 0.04062 39.6 2.171 0.0889
PS - N 0.14237 0.04062 39.6 3.505 0.0032
Results are averaged over the levels of: ROI
P value adjustment: tukey method for comparing a family of 3 estimates
Figure: time course by condition for each ROI
dat_psc_orig$ROI <- factor(dat_psc_orig$ROI, levels=c("RTPJ", "LTPJ", "PC", "DMPFC"))
dat_psc_orig$Violation <- factor(dat_psc_orig$Violation, levels=c("PS", "PH", "N"))
dat_psc_asd <- dat_psc_orig%>%
filter(Group == 'ASD' & (Violation == 'PH' | Violation == 'PS' | Violation == 'N')) %>%
group_by(Subject, Violation, ROI, Timepoint) %>%
summarise(PSC=mean(PSC))
cols <- c("PS"="darkorchid4", "PH"="red", "N"="slategray")
rois <- c(RTPJ="rTPJ", LTPJ="lTPJ", PC="precuneus", DMPFC="dmPFC")
ggplot(dat_psc_asd, aes(y=PSC, x=Timepoint, color=Violation, fill=Violation)) +
geom_smooth(na.rm=TRUE) +
facet_wrap(~ROI, ncol=4, labeller=labeller(ROI=rois)) +
annotate("rect", xmin=5, xmax=27, ymin=-Inf, ymax=Inf, alpha=.1) +
scale_x_continuous(limits=c(0,28), breaks=seq(0,28,2)) +
ylab("Percent signal change (PSC)") +
xlab("Timepoint (s)") +
scale_fill_manual(name="Condition\n", labels=c("Psychological harm", "Physical harm", "Neutral act"), values=cols) +
scale_colour_manual(name="Condition\n", labels=c("Psychological harm", "Physical harm", "Neutral act"), values=cols) +
theme_bw() +
theme(axis.text=element_text(size=16),
axis.title=element_text(size=24,face="bold"),
axis.title.y=element_text(margin=margin(r=20)),
axis.title.x=element_text(margin=margin(t=20)),
legend.text=element_text(size=20),
legend.title=element_text(size=24,face="bold"),
legend.key.size=unit(3, "lines"),
strip.text=element_text(size=28, face="bold"),
panel.grid.major.x=element_blank(),
panel.grid.major.y=element_blank())Subset data by ROI (ASD only) and define models
dat_rtpj_asd <- subset(dat_psc_long, ROI == "RTPJ" & Group == "ASD")
dat_ltpj_asd <- subset(dat_psc_long, ROI == "LTPJ" & Group == "ASD")
dat_pc_asd <- subset(dat_psc_long, ROI == "PC" & Group == "ASD")
dat_dmpfc_asd <- subset(dat_psc_long, ROI == "DMPFC" & Group == "ASD")
model_rtpj_asd <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
dat=dat_rtpj_asd, REML=FALSE)
model_ltpj_asd <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
dat=dat_ltpj_asd, REML=FALSE)
model_pc_asd <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
dat=dat_pc_asd, REML=FALSE)
model_dmpfc_asd <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
dat=dat_dmpfc_asd, REML=FALSE)rTPJ
# test pairwise contrasts
lsmeans(model_rtpj_asd, pairwise ~ Violation)$lsmeans
Violation lsmean SE df lower.CL upper.CL
PH 0.13685 0.06231 32.74 0.01003 0.2637
PS 0.15947 0.06231 32.74 0.03266 0.2863
N 0.04914 0.06231 32.74 -0.07767 0.1760
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
PH - PS -0.02263 0.0796 36.01 -0.284 0.9565
PH - N 0.08770 0.0796 36.01 1.102 0.5192
PS - N 0.11033 0.0796 36.01 1.386 0.3587
P value adjustment: tukey method for comparing a family of 3 estimates
lTPJ
# test pairwise contrasts
lsmeans(model_ltpj_asd, pairwise ~ Violation)$lsmeans
Violation lsmean SE df lower.CL upper.CL
PH 0.2884 0.08445 26.34 0.11496 0.4619
PS 0.3231 0.08445 26.34 0.14960 0.4966
N 0.1849 0.08445 26.34 0.01144 0.3584
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
PH - PS -0.03464 0.08369 420 -0.414 0.9099
PH - N 0.10352 0.08369 420 1.237 0.4321
PS - N 0.13816 0.08369 420 1.651 0.2257
P value adjustment: tukey method for comparing a family of 3 estimates
precuneus
# test pairwise contrasts
lsmeans(model_pc_asd, pairwise ~ Violation)$lsmeans
Violation lsmean SE df lower.CL upper.CL
PH 0.006022 0.0449 27.86 -0.085982 0.09803
PS 0.099761 0.0449 27.86 0.007757 0.19177
N -0.055337 0.0449 27.86 -0.147341 0.03667
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
PH - PS -0.09374 0.04262 34.72 -2.199 0.0854
PH - N 0.06136 0.04262 34.72 1.440 0.3322
PS - N 0.15510 0.04262 34.72 3.639 0.0025
P value adjustment: tukey method for comparing a family of 3 estimates
dmPFC
# test pairwise contrasts
lsmeans(model_dmpfc_asd, pairwise ~ Violation)$lsmeans
Violation lsmean SE df lower.CL upper.CL
PH 0.05501 0.06449 288 -0.071909 0.18194
PS 0.12049 0.06449 288 -0.006437 0.24741
N -0.04427 0.06449 288 -0.171196 0.08265
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
PH - PS -0.06547 0.0912 288 -0.718 0.7531
PH - N 0.09929 0.0912 288 1.089 0.5218
PS - N 0.16476 0.0912 288 1.807 0.1692
P value adjustment: tukey method for comparing a family of 3 estimates
Define the models
dat_rtpj <- subset(dat_psc_long, ROI == "RTPJ")
dat_ltpj <- subset(dat_psc_long, ROI == "LTPJ")
dat_pc <- subset(dat_psc_long, ROI == "PC")
dat_dmpfc <- subset(dat_psc_long, ROI == "DMPFC")
model_rtpj <- lmer(PSC ~ Violation*Group + (1|Subject) + (1|Item),
dat=dat_rtpj, REML=FALSE)
model_ltpj <- lmer(PSC ~ Violation*Group + (1|Subject) + (1|Item),
dat=dat_ltpj, REML=FALSE)
model_pc <- lmer(PSC ~ Violation*Group + (1|Subject) + (1|Item),
dat=dat_pc, REML=FALSE)
model_dmpfc <- lmer(PSC ~ Violation*Group + (1|Subject) + (1|Item),
dat=dat_dmpfc, REML=FALSE)
model_nt_vs_asd <- lmer(PSC ~ Violation*Group*ROI + (1|Subject) + (1|Item),
dat=dat_psc_long, REML=FALSE)Test interaction: Condition x Group x ROI
anova(update(model_nt_vs_asd, . ~ . -Violation:Group:ROI), update(model_nt_vs_asd, . ~ . -Violation:Group:ROI - Violation:Group))Data: dat_psc_long
Models:
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) +
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group): Violation:ROI + Group:ROI
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) +
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI): Violation:Group + Violation:ROI + Group:ROI
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group) 19 9420 9544 -4691 9382
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI) 21 9424 9561 -4691 9382 0.16 2 0.92
Test interaction: Condition x Group
anova(update(model_nt_vs_asd, . ~ . -Violation:Group:ROI), update(model_nt_vs_asd, . ~ . -Violation:Group:ROI - Violation:Group))Data: dat_psc_long
Models:
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) +
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group): Violation:ROI + Group:ROI
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) +
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI): Violation:Group + Violation:ROI + Group:ROI
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group) 19 9420 9544 -4691 9382
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI) 21 9424 9561 -4691 9382 0.16 2 0.92
Subset data by ROI
rTPJ
anova(model_rtpj, update(model_rtpj, . ~ . -Violation:Group))Data: dat_rtpj
Models:
update(model_rtpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
model_rtpj: PSC ~ Violation * Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_rtpj, . ~ . - Violation:Group) 7 2415 2452 -1201 2401
model_rtpj 9 2419 2466 -1201 2401 0 2 1
lTPJ
anova(model_ltpj, update(model_ltpj, . ~ . -Violation:Group))Data: dat_ltpj
Models:
update(model_ltpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
model_ltpj: PSC ~ Violation * Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_ltpj, . ~ . - Violation:Group) 7 3095 3131 -1541 3081
model_ltpj 9 3099 3145 -1541 3081 0.04 2 0.98
precuneus
anova(model_pc, update(model_pc, . ~ . -Violation:Group))Data: dat_pc
Models:
update(model_pc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
model_pc: PSC ~ Violation * Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_pc, . ~ . - Violation:Group) 7 1628 1665 -807 1614
model_pc 9 1631 1679 -807 1613 0.23 2 0.89
dmPFC
anova(model_dmpfc, update(model_dmpfc, . ~ . -Violation:Group))Data: dat_dmpfc
Models:
update(model_dmpfc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
model_dmpfc: PSC ~ Violation * Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_dmpfc, . ~ . - Violation:Group) 7 1657 1690 -822 1643
model_dmpfc 9 1661 1704 -822 1643 0.06 2 0.97
Test main effect: Condition
anova(update(model_nt_vs_asd, . ~ . -Violation:Group:ROI -Violation:Group), update(model_nt_vs_asd, . ~ . -Violation:Group:ROI -Violation:Group -Violation))Data: dat_psc_long
Models:
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) +
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group): Violation:ROI + Group:ROI
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group - Violation): PSC ~ Group + ROI + (1 | Subject) + (1 | Item) + Violation:ROI +
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group - Violation): Group:ROI
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group) 19 9420 9544 -4691 9382
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group - Violation) 19 9420 9544 -4691 9382 0 0 <2e-16
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group)
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group - Violation) ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Subset data by ROI
rTPJ
anova(update(model_rtpj, . ~ . -Violation:Group), update(model_rtpj, . ~ . -Violation:Group -Violation))Data: dat_rtpj
Models:
update(model_rtpj, . ~ . - Violation:Group - Violation): PSC ~ Group + (1 | Subject) + (1 | Item)
update(model_rtpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_rtpj, . ~ . - Violation:Group - Violation) 5 2419 2445 -1204 2409
update(model_rtpj, . ~ . - Violation:Group) 7 2415 2452 -1201 2401 7.38 2 0.025 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
lTPJ
anova(update(model_ltpj, . ~ . -Violation:Group), update(model_ltpj, . ~ . -Violation:Group -Violation))Data: dat_ltpj
Models:
update(model_ltpj, . ~ . - Violation:Group - Violation): PSC ~ Group + (1 | Subject) + (1 | Item)
update(model_ltpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_ltpj, . ~ . - Violation:Group - Violation) 5 3097 3123 -1543 3087
update(model_ltpj, . ~ . - Violation:Group) 7 3095 3131 -1541 3081 5.71 2 0.058 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
precuneus
anova(update(model_pc, . ~ . -Violation:Group), update(model_pc, . ~ . -Violation:Group -Violation))Data: dat_pc
Models:
update(model_pc, . ~ . - Violation:Group - Violation): PSC ~ Group + (1 | Subject) + (1 | Item)
update(model_pc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_pc, . ~ . - Violation:Group - Violation) 5 1636 1663 -813 1626
update(model_pc, . ~ . - Violation:Group) 7 1628 1665 -807 1614 12.5 2 0.0019 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
dmPFC
anova(update(model_dmpfc, . ~ . -Violation:Group), update(model_dmpfc, . ~ . -Violation:Group -Violation))Data: dat_dmpfc
Models:
update(model_dmpfc, . ~ . - Violation:Group - Violation): PSC ~ Group + (1 | Subject) + (1 | Item)
update(model_dmpfc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_dmpfc, . ~ . - Violation:Group - Violation) 5 1660 1684 -825 1650
update(model_dmpfc, . ~ . - Violation:Group) 7 1657 1690 -822 1643 6.95 2 0.031 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Test main effect: Group
anova(update(model_nt_vs_asd, . ~ . -Violation:Group), update(model_nt_vs_asd, . ~ . -Violation:Group -Group))Data: dat_psc_long
Models:
update(model_nt_vs_asd, . ~ . - Violation:Group): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) +
update(model_nt_vs_asd, . ~ . - Violation:Group): Violation:ROI + Group:ROI + Violation:Group:ROI
update(model_nt_vs_asd, . ~ . - Violation:Group - Group): PSC ~ Violation + ROI + (1 | Subject) + (1 | Item) + Violation:ROI +
update(model_nt_vs_asd, . ~ . - Violation:Group - Group): Group:ROI + Violation:Group:ROI
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_nt_vs_asd, . ~ . - Violation:Group) 27 9436 9612 -4691 9382
update(model_nt_vs_asd, . ~ . - Violation:Group - Group) 27 9436 9612 -4691 9382 0 0 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Subset data by ROI
rTPJ
anova(update(model_rtpj, . ~ . -Violation:Group), update(model_rtpj, . ~ . -Violation:Group -Group))Data: dat_rtpj
Models:
update(model_rtpj, . ~ . - Violation:Group - Group): PSC ~ Violation + (1 | Subject) + (1 | Item)
update(model_rtpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_rtpj, . ~ . - Violation:Group - Group) 6 2414 2446 -1201 2402
update(model_rtpj, . ~ . - Violation:Group) 7 2415 2452 -1201 2401 1.03 1 0.31
lTPJ
anova(update(model_ltpj, . ~ . -Violation:Group), update(model_ltpj, . ~ . -Violation:Group -Group))Data: dat_ltpj
Models:
update(model_ltpj, . ~ . - Violation:Group - Group): PSC ~ Violation + (1 | Subject) + (1 | Item)
update(model_ltpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_ltpj, . ~ . - Violation:Group - Group) 6 3093 3124 -1541 3081
update(model_ltpj, . ~ . - Violation:Group) 7 3095 3131 -1541 3081 0.3 1 0.58
precuneus
anova(update(model_pc, . ~ . -Violation:Group), update(model_pc, . ~ . -Violation:Group -Group))Data: dat_pc
Models:
update(model_pc, . ~ . - Violation:Group - Group): PSC ~ Violation + (1 | Subject) + (1 | Item)
update(model_pc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_pc, . ~ . - Violation:Group - Group) 6 1627 1659 -807 1615
update(model_pc, . ~ . - Violation:Group) 7 1628 1665 -807 1614 1.16 1 0.28
dmPFC
anova(update(model_dmpfc, . ~ . -Violation:Group), update(model_dmpfc, . ~ . -Violation:Group -Group))Data: dat_dmpfc
Models:
update(model_dmpfc, . ~ . - Violation:Group - Group): PSC ~ Violation + (1 | Subject) + (1 | Item)
update(model_dmpfc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_dmpfc, . ~ . - Violation:Group - Group) 6 1655 1684 -822 1643
update(model_dmpfc, . ~ . - Violation:Group) 7 1657 1690 -822 1643 0 1 0.95
Data organization
# transform data to long format
dat_mvpa_long <- dat_mvpa_orig %>%
gather(cond, acc, which(regexpr("_", colnames(.)) > 0)) %>%
separate(cond, c("roi", "type"), extra="drop")Comparing mean accuracies to chance (50%)
mvpa_res <- matrix(nrow=0, ncol=6,
dimnames=list(NULL, c("group","col", "t", "df", "p.value", "estimate")))
groups <- levels(dat_mvpa_orig$group)
for (cl in 3:ncol(dat_mvpa_orig)) {
for (grp in groups) {
temp <- dat_mvpa_orig[dat_mvpa_orig$group == grp, cl]
mvpa_res <- rbind(mvpa_res, c(grp, colnames(dat_mvpa_orig)[cl],
unlist(t.test(temp, mu=.5, alternative="greater"))[c(1, 2, 3, 6)]))
}
}
# tidy up results
mvpa_res <- tidy(mvpa_res)
mvpa_res[,3:6] <- lapply(mvpa_res[,3:6], function(x) as.numeric(as.character(x)))
mvpa_resComparing mean accuracies of experimental tests vs. permutation tests
mvpa_res_exp_vs_perm <- matrix(nrow=0, ncol=8,
dimnames=list(NULL, c("group","col1", "col2", "t", "df", "p.value", "estimate.of.exp","estimate.of.perm")))
groups <- levels(dat_mvpa_orig$group)
for (cl in seq(3, ncol(dat_mvpa_orig), by=2)) {
for (grp in groups) {
temp1 <- dat_mvpa_orig[dat_mvpa_orig$group == grp, cl]
temp2 <- dat_mvpa_orig[dat_mvpa_orig$group == grp, cl+1]
mvpa_res_exp_vs_perm <- rbind(mvpa_res_exp_vs_perm, c(grp, colnames(dat_mvpa_orig)[cl], colnames(dat_mvpa_orig)[cl + 1],
unlist(t.test(temp1, temp2, alternative="greater"))[c(1, 2, 3, 6, 7)]))
}
}
# tidy up results
mvpa_res_exp_vs_perm <- tidy(mvpa_res_exp_vs_perm)
mvpa_res_exp_vs_perm[,4:8] <- lapply(mvpa_res_exp_vs_perm[,4:8], function(x) as.numeric(as.character(x)))
mvpa_res_exp_vs_permComparing mean accuracies of NT group vs. ASD group
mvpa_res_nt_vs_asd <- matrix(nrow=0, ncol=6,
dimnames=list(NULL, c("col", "t", "df", "p.value", "estimate of NT", "estimate of ASD")))
groups <- levels(dat_mvpa_orig$group)
for (cl in 3:ncol(dat_mvpa_orig)) {
temp1 <- dat_mvpa_orig[dat_mvpa_orig$group == 'NT', cl]
temp2 <- dat_mvpa_orig[dat_mvpa_orig$group == 'ASD', cl]
mvpa_res_nt_vs_asd <- rbind(mvpa_res_nt_vs_asd, c(colnames(dat_mvpa_orig)[cl],
unlist(t.test(temp1, temp2, alternative="two.sided"))[c(1, 2, 3, 6, 7)]))
}
# tidy up results
mvpa_res_nt_vs_asd <- tidy(mvpa_res_nt_vs_asd)
mvpa_res_nt_vs_asd[,2:6] <- lapply(mvpa_res_nt_vs_asd[,2:6], function(x) as.numeric(as.character(x)))
mvpa_res_nt_vs_asdFigure: classification accuracies
dat_mvpa_plot <- dat_mvpa_long
dat_mvpa_plot$roi <- factor(dat_mvpa_plot$roi, levels=c("RTPJ", "LTPJ", "PC", "DMPFC"))
dat_mvpa_plot$group <- factor(dat_mvpa_plot$group, levels=c("NT", "ASD"))
ggplot(dat_mvpa_plot, aes(y=acc, x=type)) +
stat_summary(fun.data="mean_cl_boot", position=position_dodge(0.4)) +
facet_grid(group ~ roi, labeller=labeller(roi=c(RTPJ="rTPJ", LTPJ="lTPJ", PC="precuneus", DMPFC="dmPFC"))) +
geom_hline(yintercept=.50) +
ylab("Classification accuracy") +
scale_x_discrete("Test type\n", labels=c("exp"="Experimental", "perm"="Permutation")) +
theme_bw() +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=16,face="bold"),
axis.title.y=element_text(margin=margin(r=20)),
axis.title.x=element_text(margin=margin(t=20)),
strip.text=element_text(size=16, face="bold"),
panel.grid.major.x=element_blank(),
panel.grid.major.y=element_blank())